#include "general.h"
#include <stdio.h>
#include <string>

void skip_line( FILE * fp );
double next( FILE * fp );
void nrerror( char * str);
void rest();
void noMoreMemory();

double *vector(int nl, int nh);
int *ivector(int nl, int nh);
double **matrix(int nrl, int nrh, int ncl, int nch);
double ***tensor(int nrl, int nrh, int ncl, int nch, int ndl, int ndh);
void free_vector(double *v, int nl, int nr);
void free_ivector(int *v, int nl, int nr);
void free_matrix(double **m, int nrl, int nrh, int ncl, int nch);
void free_tensor(double ***t, int nrl, int nrh, int ncl, int nch, int ndl, int ndh);
void bandec(double **a, int n, int m1, int m2, double **al, int indx[], double *d);
void banbks(double **a, int n, int m1, int m2, double **al, int indx[], double b[]);
void tridag(double a[], double b[], double c[], double r[], double u[], unsigned long n);
void lubksb(double **a, int n, int *indx, double b[]);
void ludcmp(double **a, int n, int *indx, double *d);
void sor(double **a, double **b, double **c, double **d, double **e, double **f, double **u, int jmax, double rjac);
void gaussj(double **a, int n, double **b, int m);

double Bernoulli( double x );
double chem_composition(double,double);
void analytical(double V, double x0, double xw, double mass0, double massqw, char * file1, char * file2);
double * fluctuation(int type, int num, double mean, double perc );


inline double FD( double x ) { // Fermi Integral function
  return 1.0/(exp(-x)+0.75*1.7725*pow((x*x*x*x)+50+33.6*x*(1.0-0.68*exp(-0.17*(x+1.0)*(x+1.0))),0.375));
}
inline double dFD( double x ) { // Derivative of the Fermi integral function
  return  1.05*pow(FD(x),2)*(exp(-x)-1.138*0.375*0.75*1.7725
			     *pow((x*x*x*x)+50+33.6*x*(1.0-0.68*exp(-0.17*(x+1.0)*(x+1.0))),-0.625)
			     *3.0*(x*x*x)+33.6*(1.0-0.68*exp(-0.17*(x+1.0)*(x+1.0)))
			     +33.6*x*0.68*0.17*2.0*(x+1.0)*exp(-0.17*(x+1.0)*(x+1.0)));
}
inline double GX( double x, double g ) { return (1.0/(1.0+(1.0/g)*exp(x)));}
inline double dGX( double x, double g ) { return  -(1.0/g)*exp(x)/pow( 1.0 + (1.0/g)*exp(x),2 );}
inline double GD( double x ) { return (1.0-(1.0/(1.0+0.5*exp(-x))));}
inline double GA( double x ) { return 1.0/( 1.0 + 4.0*exp(-x) );}
inline double dGD( double x ) { return - (0.5*exp(-x))/pow(1.0+0.5*exp(-x),2);}
inline double dGA( double x ) { return  4.0*exp(-x)/pow( 1.0 + 4.0*exp(-x),2 );}

